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We study a model of a composite system constructed from a "pairing layer" of disconnected 
attractive-(7 Hubbard sites that is coupled by single-particle tunneling, t±, to a disordered metallic 
layer. For small inter-layer tunneling the system is described by an effective long-range XY phase 
model whose critical temperature, T c , is essentially insensitive to the disorder and is exponentially 
suppressed by quantum fluctuations. T c reaches a maximum for intermediate values of t±, which 
we calculate using a combination of mean-field, classical and quantum Monte Carlo methods. The 
maximal T c scales as a fraction of the zero temperature gap of the attractive sites when U is smaller 
than the metallic bandwidth, and is bounded by the maximal T c of the two-dimensional attractive 
Hubbard model for large U. Our results indicate that a thin, rather than a thick, metallic coating 
is better suited for the enhancement of T c at the surface of a phase fluctuating superconductor. 
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I. INTRODUCTION 

Recently, composite systems made of a metal over- 
laying a superconductor with low phase stiffness have 
received theoretical attention 1 - - ™ after experiments have 
demonstrated that the superconducting transition tem- 
perature, T c , can be enhanced in bilayers consisting of 
underdoped and overdoped cuprates. 4 ' 5 In particular, 
Ref. [l| considered a model of a pairing layer with zero 
phase stiffness, constructed from disconnected attractive- 
U Hubbard sites, that is coupled via single-particle tun- 
neling t± to a metallic noninteracting layer. As tj_ in- 
creases from zero, phase coupling between the pairing 
sites is established by Josephson tunneling through the 
metallic layer. At the same time, however, the pairing 
strength is diminished owing to the proximity effect in- 
duced by the same derealization events. Consequently, 
T c reaches a maximum at an intermediate value of t±, 
where both these trends attain a simultaneous optimum. 
It was shown, within a mean-field approximation, that 
the maximal T c approaches the limit set by the pairing 
scale of the Hubbard sites when the metallic bandwidth, 
Ht, is much larger than U. Here we wish to reexamine 
these results using more accurate analytical and numer- 
ical techniques. 

An additional goal of the present work is to study the 
effects of disorder in the metallic layer on T c . Owing to 
Anderson- , it is well known that T c of a dirty BCS s- 
wave superconductor is essentially insensitive to disorder 
as long as T c is much larger than the local mean level 
spacing. However, BCS theory ignores fluctuations in 
the superconducting order parameter and one expects 
deviations from Anderson's result when fluctuations are 
large. Indeed, phase fluctuations dominate the physics 
in strongly disordered and granular systems where 
they can reduce T c to zero. In the bilayer that we study 
the disorder resides away from the pairing sites and is 
therefore expected to have only mild consequences for 
pairing. However, since phase stiffness is established by 



Josephson tunneling through the metal it is interesting 
to study the manner in which it is affected by the random 
potential. 

Our results indicate that although the general picture 
obtained in Ref. [l| holds, it requires some modifica- 
tions. First, the analysis of Ref. [lj relied on a Bogoli- 
ubov de-Gennes (BdG) mean-field treatment to calcu- 
late the phase stiffness. While this approximation takes 
into account the suppression of the phase stiffness by 
thermally excited quasiparticles, it ignores the renormal- 
ization of the stiffness by thermal and quantum phase 
fluctuations. The latter have little effect in the nearest- 
neighbor XY model and the two-dimensional attractive 
Hubbard model^, where using the unrenormailzed stiff- 
ness reproduces the correct T c to within 40%. We show 
that this is not the case in the presence of long range 
phase couplings. Such couplings, which extend up to the 
thermal length of the metal, occur in the small t ± regime 
of the bilayer. Under these circumstances classical phase 
fluctuations on scales smaller than the thermal length 
lead to a rapid decrease of the stiffness and therefore to 
a much lower T c than is anticipated from the unrenor- 
malized stiffness. Even more important are the quantum 
phase fluctuations in the small t± regime, which lead to 
an exponential suppression of T c relative to its mean- field 
value. 

Secondly, using quantum Monte Carlo and classical 
Monte Carlo mean-field techniques we find that the high- 
est T c (maximized over t±), is smaller by a factor 3-4 than 
the mean- field prediction for small and intermediate U/t. 
For large U jt it is bounded by the maximal T c of the two- 
dimensional attractive Hubbard model. Our calculations 
reveal that the maximum is largely governed by classical 
phase fluctuations. 

For small inter-layer tunneling the primary effect of 
the disorder is to decrease the range over which coherent 
phase coupling is mediated through the metal. However, 
as long as the metal remains in the diffusive regime, this 
has a negligible effect on T c . Using a mean- field approach 
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to estimate the effects of disorder away from the small 
t± regime yields a weak correction to the maximal T c at 
small disorder strength. Nevertheless, once the disorder 
becomes large in comparison to the hopping amplitude 
in the metallic layer the maximal T c is significantly sup- 
pressed. We conclude the paper with a short discussion 
of the relevance of these insights to attempts to enhance 
T c at the surface of a phase fluctuating superconductor. 



II. THE MODEL AND ITS ANALYSIS IN THE 
SMALL t± LIMIT 

A. The model 

We consider a bilayer, see Fig. [T] consisting of a nonin- 
teracting disordered upper layer and a lower layer of dis- 
connected negative-?/ Hubbard sites. Each layer contains 
N = L 2 /a 2 sites in a square array with lattice constant 
a. Neighboring sites on the two layers are connected via 
single particle tunneling, t±. We denote by Ci„ and fi„ 
the annihilation operator of an electron with spin a on 
the ith site of the upper and lower layer, respectively. 
The imaginary time action describing the model is 
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Here /3 — l/T is the inverse temperature, [i is the chem- 
ical potential and t is the hopping amplitude between 
sites in the upper layer, which contains a Gaussian ran- 
dom potential V and a constant potential e used to adjust 
the Fermi energy away from any van-Hove singularities. 



B. Phase-only action 

We proceed by decoupling the interaction term using a 
Hubbard- Stratonovich transformation —Uf^f^fufa — > 

AfAi/tf + Alfijit + Ajyi, where A, = |A,.| e *«. 
Gauge transforming f ia — > fi a e l6i ^ 2 , integrating out the 
noninteracting layer and denoting = leads 
to 

S = [ drdr' V *t(T)[M<°> + ^ (0) + F (1) ]^(r') 

+ [ drJ2\^(r)\ 2 , (2) 
Jo 




FIG. 1: The model: A square lattice of disconnected 
attractive-C/ Hubbard sites is coupled via single-particle tun- 
neling t± to a noninteracting tight-binding layer with on-site 
random disorder. 



with 



= [d T I - ficr 3 + lA^lo-i^^T - r'), 



-9 T 6»i(r)cr3%(5(T - t'), 

W-^(-')1/2^.. (t _ t / )<73) (3) 



t 2 , e 



where Qij (r) is the Green's function of the noninteracting 
layer, and er are the Pauli matrices. 

Next, we integrate out the degrees of freedom of the 
pairing layer. We do so perturbatively to second or- 
der in V~(°' and The time independent amplitude 
| Ai| = |Aj|(O n = 0) is essentially set by the zeroth 
order contribution {0/U)J2i l A *| 2 - TrlnM(°)(|Ai|) to 
the effective action, which provides within a saddle point 
approximation the BCS gap equation for the decoupled 
Hubbard sites. We assume that the pairing layer is close 
to half filling such that fi <C |A|, and obtain as a re- 
sult the T = solution |A| = A ~ U/2. The 0(t\) 
contribution modifies the gap equation and leads to a 
reduction of order t\/t of Ao, reflecting the proximity 
effect. In the following we assume that t\Jt,T <C Ao 
and therefore neglect both this effect and any thermal or 
quantum fluctuations of the gap amplitude. 

The action governing the phase fluctuations is derived 
in the appendix. There we show that its dominant terms 



Se = 



1 



8Ao j,, 



dTY J {dr9 l f 



f drdr' V K( nj , r - r') cos[^ (r) - Bj (r')] . (4) 
Jo 



The Kernel K(r,r) decays exponentially for spatial sep- 
arations Tij = |rj — Tj\ larger than the thermal length 
It- This decay reflects the loss of coherence between the 
dynamical phases of the members of a pair that medi- 
ates the phase coupling. Due to thermal smearing of the 
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Fermi-Dirac distribution the difference between the en- 
ergies of the two electrons is of order T, which leads to 
loss of coherence after time \/T . In the clean limit where 
the elastic mean free time of the metal satisfies t^ 1 <C T 
this coherence time translates to a distance It = vp/T 
covered by the ballistically propagating electrons with 
Fermi velocity v F . In the same clean limit we find for 
It T » £ = vp/A , that the kernel is 

, n t 4 ,N F a 4 l 1 



2tt% f r (A r) 2 + (r/0 2 

where TVf is the density of states of the metallic layer 
at the Fermi energy. Up to logarithmic corrections the 
behavior at shorter distances may be approximated by 

- 4 1 



K(£ > r > a, r) = 



!L p -2A |r| 



(Aor) 



1 + (Aor) 4 



(6) 



Owing to the reasons outlined above the phase cou- 
pling in the diffusive regime r" 1 3> T exhibits a similar 
exponential decay beyond the thermal length, which for 
a disordered metal with diffusion constant D is given by 
It = \/D/T. In the important range It ^ r ^> £ = 
^/D/ Aq the kernel behaves to within logarithmic accu- 
racy as 



K(l T > r > £,t) 



d JV F a 4 



27r 2 ^ (A r) 2 + (r/2£) 4 ' 



(7) 



with phase couplings /(r) given by 



£ » r > a, Z e 
Zt > r > £ 
r > Z T 



clean limit 

27T V r/ 



2a 2 -2irr/h 
ItT 



diffusive limit 

Kf) 2 '°(^f) 



1 f 

7r V r y 



(9) 



/_a_\ ' (27T) 1 /' 



where Z e = UfT e is the elastic mean free path and 7 is 
the Euler constant. In our region of interest T <C Ao 
and the physics is governed by the couplings in the range 
It ^> T ^> £■ Thus, we are led to study the following 
XY-type model 



e~ n ' /x cosf 



9 3 ), (10) 



with J ~ t\Npa 2 1 A\ and A ~ It as an effective descrip- 
tion of the classical phase fluctuations. Here, for sake of 
simplicity, we extended the r~ 2 behavior of the coupling 
down to the short distance cutoff, ignoring the crossover 
when r < £, see Eq. ©. As we shall demonstrate, this 
leads to a negligible correction to T c . 

Since the phase ordering transition belongs to the BKT 
class its critical temperature is related to the universal 
jump of the renormalized phase stiffness at criticality: 



C. Classical phase fluctuations 

In the small t± limit the superconducting critical tem- 
perature, T c , equals the phase ordering temperature as 
determined by the action, Eq. ((4]). Being a finite tem- 
perature transition in a two-dimensional system with fi- 
nite range couplings it is clear that the phase ordering 
transition belongs to the classical Berezinskii-Kosterlitz- 
Thouless (BKT) universality class. However, T c itself is 
determined by both quantum (time-dependent) and clas- 
sical (time-independent) phase fluctuations. As we shall 
see, the quantum phase fluctuations can not be neglected 
and are in fact dominant for small t±, leading to expo- 
nential suppression of T c . Albeit, we begin by considering 
the effects of the classical fluctuations. We do so since 
the treatment of the classical fluctuations will reveal a 
lesson concerning the renormalization of the phase stiff- 
ness which is generic to models with interactions that 
extend well beyond nearest neighbors and since it will 
provide the tool to calculate the effects of the quantum 
fluctuations. 

To this end, consider the case of a time-independent 
but space fluctuating 9. Assuming /3Aq 3> 1 and carrying 
out the time integration in Eq. (j4]) results in 



Ps{T c ) 



(11) 



-0 



t\N F a 2 

A 2 
^0 



W, (8) 



In turn, the phase stiffness is calculated from the free 
energy in the presence of a phase twist, </>, per bond in 
the x direction. That is, if H = —(1/2) J^i j Jij cos[#i — 
8j + (xi — Xj)4>/a] then 



Ps{T) 




Xj) 2 Jij cos( 



x j) Jij sin( 



(12) 



where here () denotes thermal averaging. 

For the standard XY model it is known that the BKT 
criterion (|11[) yields a fair estimate for T c even when p s 
is replaced by the bare stiffness unrenormalized by 
vortices and longitudinal phase fluctuations. Indeed, for 
that model p° s = J, as calculated from Eq. (|12p using a 
uniform phase field 9i = 6. This gives the estimate T c = 
(n/2)J, which is to be compared with the most recent 
numerical value 15 T c = 0.8929J. One may, therefore, 
attempt to apply the same approximation to calculate 
the transition temperature of model ([IT))) . For this case 
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FIG. 2: The phase stiffness of the effective model, Eq. ([TDJ, 
with A/a = 4. The dashed line is a fit (with A/a = 3.897, 
J/J = 0.108 and 2 = 8) using the mean-field approximation 
discussed in the text. The inset shows the size dependence 
of p a in the vicinity of the transition, which occurs near the 
crossing of the stiffness curves with the solid black line p s = 
(2/n)T. 



one finds, using the fact A»a, that 

rf-s/*'®'^ -S (:)"'■ < 13 > 

and the estimate T c = (7rA/2a) 2 J. Accordingly, the tem- 
perature dependence of A ~ It would imply T c ~ i^/Ao 
for the diffusive system and T c ~ (t^Npvp/ A 2 ,) 1 / 3 in the 
clean limit. The latter result was previously derived in 
Ref. [l| within the same approximation. 

To test the validity of these estimates we calculated 
p s {T) for the Hamiltonian, Eq. (ITUl) . via Monte-Carlo 
simulations of systems with up to L = 512. For this 
purpose we used Eq. ([12)) and implemented the Wolff 
algorithm^, which is easily generalized to include cou- 
plings beyond nearest-neighbor range. Our results, see 
Fig. [3J indicate that the phase stiffness develops a discon- 
tinuity, in accord with the expected signature at a BKT 
transition. However, unlike the situation in the standard 
XY model, p s (T c ) is massively renormalized down from 
its bare (T — 0) value p° s . In fact, as shown by Fig. [3l this 
renormalization leads in our region of interest A/(i> 1 to 
T c ~ 3.21Jln(0.74A/a) instead of T c ~ J(A/a) 2 . There- 
fore, using the above stated values of J and A in terms of 
the parameters of the original model we find the following 
estimate for T C) based on thermal fluctuations alone 

N F a 2 tj ^ VTV^Aq A, (A 
T c , class ~ In )+0(t x lnlni ± ) , 

,( 14 ) 

where 77 = D,vpa in the diffusive and ballistic regimes, 
respectively. 




1 fO 20 

A/a 



FIG. 3: The critical temperature of the effective model, Eq. 
(|10p. as function of A/a. 

In order to acquire an insight into these findings let us 
consider the following coarse-grained model where we di- 
vide the system into blocks containing (A/a) 2 unit spins. 
Each spin is assumed to interact with all the other spins 
in its own block and in the z neighboring blocks. To 
make analytical progress we replace the original coupling 
in our model, which decays as 1/r 2 , by its average over 
a block. Hence, if we denote by s|, i = 1 . . . (A/a) 2 , the 
spins on the 7th block, and by mj = Yli s f t ne total 
(super)spin on that block, we are led to study the model 

H = -^m 2 - j^mrmj, (15) 

1 (I, J) 

where (I, J) denotes the blocks that couple to block I 
and 




The coarse-grained model, Eq. (TT5|) . is a model of soft 
superspins whose average length M — (|m|) is deter- 
mined by the intra-block phase fluctuations. As long as 
fluctuations in M are ignored the system can be viewed 
as an ordinary XY model for unit superspins with cou- 
pling JM 2 . Hence, we are interested in calculating the 
temperature dependence of the latter. To this end, we 
continue to treat the intra-block fluctuations exactly but 
treat the inter-block coupling in mean-field approxima- 
tion using the following effective single block Hamiltonian 

H MF = ~'^m. 2 - JzM-m, (17) 

with the self-consistency condition M = (rn)g- . 
Using the Hubbard-Stratonovich transformation^ 

i r°° — - 

e 0J m */2 = i / d 2 ye ^ + ^j m -y^ (18) 
J -00 
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we are able to write the mean-field partition function as 

-y 2 + |x|cos0^ Qgs 



J MF 



oo >2 ( A / a ) /-27T 

n 



n 



i=i 



where cf>i is the angle between Si and x = y/20Jy + 
z/3JM. Carrying out the integrations over the </>;S and 
the direction of x we arrive at 



'MF 



(2tt 



,(A/a) 2 



(3J 

poo 2 

x / dxx[I (x)] (X/a) I Q (zMx)e- x2/(2fsJ \ (20) 
Jo 

where Iq(x) is the modified Bessel function of the first 
kind. Finally, the self-consistency condition implies 
M 2 = (f3J)~ 1 d\nZMF/dz, from which follows 



M 



1 J^dxx 2 [I (x)} {X/a) h(zMx)e~ x2 /W J ^ 
(l + z)pj J™dxx[I (x)] iX/a)2 I (zMx)e-* 2 /Wj) ' 

Near the mean- field transition temperature, T c- mf, 
M — > and we may approximate Iq(zMx) ~ 1 and 
Ii(zMx) « zMx. This, together with the parameteriza- 
tion /3 C J = c(a/A) 2 for the inverse critical temperature, 
turns Eq. ((21} into 



(A/a) 2 



1 



z Jo 00 ^ 3 [l {V2~c{a/\)x)} 



(22) 



1 + 2 r °° [I (y/2c{a/X)x)] ' 

Since 1q(:e -C 1) ~ 1 + x 2 /4 one finds in the limit a/A <C 

1 that [l (V2^(a/X)x)] {X/a)2 w e ca;2 / 2 . The integrals 
in Eq. (1221) are then readily evaluated with the result 
c = 2/(1 + z), leading to 



l + z /A 



,MF 



J ~ Jln(A/a). 



(23) 



The phase stiffness of the coarse-grained model is de- 
termined by two types of processes: intra-block fluctu- 
ations that reduce M and with it the average coupling 
between superspins, and inter-block fluctuations of the 
superspins. Our mean-field treatment ignores the second 
type of fluctuations. Fig. [2] depicts a fit to p s (T)/J of 
Hamiltonian (TTU)) using M 2 (T) obtained by solving Eq. 
([21]) . The fit begins to deviate from p s {T) as the latter 
becomes of order (2/tt)T. Hence, the following physical 
picture emerges: The phase stiffness is rapidly reduced 
from its large bare value by fluctuations on scales smaller 
than A. These include both longitudinal and transverse 
vortex excitations which reach, according to our numer- 
ical findings, much higher densities below T c as com- 
pared to the case of the standard XY model. It is this 
rcnormalization that is responsible for the T c scaling with 
ln(A/a). Once p s approaches (2/w)T vortex fluctuations 
on scales larger that A become important and drive the 
system through a BKT transition. 



Before proceeding to discuss the quantum fluctuations 
let us consider the approximation we made in neglecting 
the crossover to a slower decay of the coupling, Eq. ^ , 
for r < £. Based on the insights gathered above we can 
easily modify the mean-field treatment to include this 
crossover by defining the average coupling according to 



d 2 r 



,-r/X 



-r/\ 



"J 



(24) 



Since A 3> £ this introduces only a small correction to 
the mean-field transition temperature, Eq. (|23p. 



D. Quantum phase fluctuations 

For small t± the short time phase dynamics on a single 
site is governed by the first term in the action (j4]) , which 
implies 



-i[9i(T)-6i(T')]\ 



-2A |t-t'| 



(25) 



This means that a site phase is essentially constant over 
a time of the order of Aq 1 and allows us to coarse grain 
space-time into "needles" of length A^" 1 in the imagi- 
nary time direction. The phases of the needles interact 
according to the last term in Eq. Q , resulting in a coarse 
grained phase action 

= g(rn , r - r') cos[0 4 (r) - 6, (r')] , (26) 



where g(r,r) = Aq 2 K(t, t), and Oiij) denotes the phase 
of the needle centered around time r at site i. 

The fact that g{r, r) is long ranged in the time direction 
allows us to apply a similar mean-field approach to the 
one employed above in order to estimate T c . Now, how- 
ever, we include the effects of both quantum and thermal 
fluctuations. For this purpose we divide space-time into 
rods of length /3 in the time direction and spatial area 
l\. Each phase within a rod is taken to interact with all 
the other phases in its rod and in the z neighboring rods 
with a coupling strength which is the average of g(r, r) 
over a rod 



13/2 



rod a /n ^, 



1 2nt 4 A _N F 

iVrod A 3 , 



dr rf(r) 



(27) 



Here A^ ro d = /3Ao(Zt/o) 2 is the number of needles within 
a rod and we used the fact that the averaging along the 
time direction leads to the couplings f(r), Eq. ©, en- 
countered previously in the context of the thermal fluc- 
tuations. According to the mean-field analysis preceding 



G 



Eq. (l23t criticality occurs when N ro dg = 2/(1 + z). As a 
result one obtains 



T c — min(A , — ) exp 
a 



2Ag 



(1 + z)t\N F a 2 



(28) 



both in the clean and diffusive limits. We conclude that 
for t± <C y/tAg, where our treatment applies, quantum 
phase fluctuations induce an exponential suppression of 
T c from its value based on thermal fluctuations only, Eq. 
(1141) . Secondly, within our approximation the disorder 
has no effect on T c as long as the metallic layer is in the 
diffusive regime. 



III. LARGE AND INTERMEDIATE t ± 
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A. The large t± limit 

When t± U,t the physics is dominated by the inter- 
layer tunneling and the energy is minimized by the cre- 
ation of a symmetric state on each c — / dimer (assuming 
that the system is below half filling). Denoting by ai a 
the annihilation operator of this state on site i one finds 
Cia ~ fia ~ a-ia/V^ and a disordered Hubbard model as 
the effective Hamiltonian 

H = -\ (4^ + H.c) 

s^fe + Vi \ U v-^ 

+ 2^ I — M - *x ) n ia -JZ^ n it n ih ( 29 ) 

icr ' i 

where ni a — a\ a ai a . In the weak coupling limit U <C t 
the phase stiffness is largely determined by the ampli- 
tude of the order parameter. Hence, the BKT critical 
temperature is very close to the BCS mean field transi- 
tion temperature^ and Anderson's theorem applies. For 
stronger interaction U ~ t the model was investigated 
using both the mean-field approximation™ and quan- 
tum Monte-Carlo (QMC) simulations^. These studies 
demonstrated how with increasing disorder strength the 
system becomes dominated by phase fluctuations which 
eventually turn it into an insulator. In the limit U 3> t 
the model can be mappe d 20 ' 21 to a nearest-neighbor 
quantum XY model whose T c ~ t 2 /U. 



B. The intermediate t± regime 

The preceding analysis shows that T c rises with small 
t± according to Eq. (|28[) . while approaching an asymp- 
totic large t± value which scales as t exp(— 16t/U) and 
t 2 /U in the weak and strong interacting limits, respec- 
tively. The absence of a small parameter in the interme- 
diate regime t± ~ U makes it a more difficult theoretical 
challenge and one needs to resort to numerical methods. 
For this purpose we took advantage of the fact that the 



FIG. 4: The critical temperature as function of t± for a clean 
bilayer with U/t = 1 and n = 1.5, calculated using QMC and 
MCMF methods. The inset depicts a comparison of these 
results with the predications of the BdG mean-field theory. 

bilayer model is free from the sign problem at all dop- 
ing levels, and implemented a determinantal QMC tech- 
nique to calculate its phase stiffness. Consequently, T c 
was evaluated via the BKT criterion. The phase stiff- 
ness was extracted from imaginary time current-current 
correlations according to a theorem by Scalapino, White 
and Zhang 2 ^ 

Ps = -\ i( K x) + & xx {q x =0,q y ^ 0,iu n = 0)] . (30) 
Here 

Kx = H (/'.•• + 4, CT c r+i,a) , (31) 

r,<7 

is half the kinetic energy per site, and in practice the 
limit q y — > of the correlation function 

A xx (q,iu, n ) = - dTe w ^(j x (ci,T)j x (-q,0)), (32) 
where 

j x (q) =itJ2 e- iq ' r (c\ + . ^ a - 4,^+^) , (33) 

r , ct 

stands for its value for q y — 2n/L in the finite systems 
that we simulate. We used the BSS algorithm^ to carry 
out the evolution in configuration space and the Hirsch 
method 2 ^ for the measurement of A xx . We found that in 
the range U > t it was sufficient to set the QMC time 
slice to At = (4J7) -1 and sweep through the system 
5000-10,000 times in order to limit the systematic and 
statistical errors of the calculated T c to a few percents. 

Our results for a clean bilayer with U/t = 1, total 
density n = 1.5 and L/a = 12, 14 appear in Fig. @] A 



7 



.12 



- 


.08 


a 










.04 









j 






H 


2 


R 






1 









1.5 




1 


< 


.5 










T ft 

c, 2D, max 



t A /2t 



< 



O MCMF 
O QMC 



H 1 1 1 1 1 1 1 H 



H 1 1 1 1 1 1 1 1 h 



XX- 



-o- 



._()_. 



2D, max 



8 12 
U/t 



16 20 



FIG. 5: Upper panel: The maximal critical temperature as 
function of U/t for a clean 12 x 12 bilayer with n = 1.5. 
The horizontal dashed line depicts the highest T c (maximized 
over U/t) of the single-layer attractive Hubbard model with 
n = 0.5. The left dashed line corresponds to the mean-field 
transition temperature of the disconnected Hubbard sites. 



tion of U/t. The dashed line is a fit to tx jmax oc ytU. Lower 
panel: A ef r(tx,max) = £x,max/Ao as function of U/t. The 
dashed line corresponds to the average amplitude of the or- 
der parameter at the highest T c of the single-layer attractive 
Hubbard model with n = 0.5. 



maximum in T c as function of t± is observed at ij_.max — 
0.6t. While this behavior is akin to the BdG mean-field 
findings of Ref. [l], shown in the inset, the two sets of 
data differ quantitatively. The peak in T c , as calculated 
by QMC, occurs at somewhat higher values of t± and 
its magnitude, T C!inax , is about 3 times smaller than the 
corresponding maximum in the BdG mean-field T c . 

The use of QMC to study the U/t dependence of T C:WLax 
is restricted by the low temperatures that are encoun- 
tered in the small U/t regime and the short At required 
when U/t is large. To partially overcome these difficul- 
ties we employ an approximate method which neglects 
quantum fluctuations of the order parameter. As we saw 
in Sec. Ill CI such an approximation fails for small t±. 
However, our findings demonstrate that it is sufficient 
near the maximal T c . In this Monte Carlo Mean Field 
method (MCMF) , (originally introduced by Mayr et al^- 



in its d-w&ve version), a classical Monte-Carlo scheme is 
used to average over random configurations of the local 
pairing field. The partition function is given by 



N 



Z = ff[dAidA*e-^^^^ Ai \ 2 Z({Ai}), (34) 

where Z({Ai}) is the partition function of the BdG quasi- 
particles for a given configuration of the pairing field. 
Observables, such as the phase stiffness, are calculated 
using the quasi-particle Green's functions and averaged 
over the pairing field configurations. The critical temper- 
ature is determined using Eq. (|3U)l . just as in the QMC 
simulations. As shown in Figs. 0] and [5] MCMF repro- 
duces to within 20% the QMC results for T Cimax over the 
range \<U/t<% (we attribute the fact that the MCMF 
T c curve lies below the corresponding QMC curve in the 
U/t = 1 system to a different finite size scaling of the two 
methods) . We therefore conclude that the maximal T c is 
largely governed by thermal fluctuations. 

Utilizing the MCMF method to calculate T Ci , nax be- 
yond the U/t range which is amenable to QMC simula- 
tions one obtains the following behavior, depicted in Fig. 
[5] First, at small U /t the maximal critical temperature 
scales with the mean- field temperature, Tmf ~ U/4, of 
the disconnected Hubbard sited. The latter is the tem- 
perature at which the pairing gap on each site closes and 
thus sets a maximum conceivable value for T c which takes 
full advantage of the pairing scale. Numerically, we find 
for U < t that T Cymax w Tmf/3. This is to be contrasted 
with the mean-field result of Ref. Q], T c max — > Tmf in 
the same range of parameters. Presumably, the differ- 
ence stems from the absence of (predominantly classical) 
phase fluctuations in the mean-field calculation. 

At large U /t, T Cimax seems to saturate to a limiting 
value, which is close to the maximal T c of the quarter 
filled two-dimensional attractive Hubbard model^ Our 
MCMF results indicate that amplitude fluctuations of the 
order parameter play a minor role at the T c maximum of 
both models, changing it by about 5%. We would there- 
fore attempt to understand the above behavior from the 
perspective of thermal phase fluctuations only. To this 
end we decouple the interaction term in the original ac- 
tion ([T]) and integrate out the pairing sites while allowing 
only spatial phase fluctuation in the pairing field, i.e. as- 
suming Aj(r) = Aoe l0i . The result is 



x c ja (u n ) + fj A off E 



t]_(iio n - n)8i 



\ 2 



e iei ct t (o;„)ct | (-a;„)+H.c. 



where 



A 



off 



f 2 
__L 

A ' 



(35) 



(36) 
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As shown in Fig. [5j ij_ !max scales as \ftU for U > t. 
Since A ~ U/2 this implies that at t± ymax and for U 3> t 
the 0(t\) correction to the first term in the action (|35|) 
can be neglected in comparison to . Moreover, since 
the important contribution to the action comes from fre- 
quencies within the metallic bandwidth the prefactor in- 
side the second term of Eq. (|55f can be replaced in the 
limit [/ > t by 1. Within these approximations Eq. 
(|35|) becomes the action of a two-dimensional attractive 
Hubbard whose filling is set by the filling of the metal- 
lic layer and whose pairing field amplitude is given by 
A c ff. The lower panel of Fig. [5] demonstrates that in 
the U > t regime A c g(t± . max ) coincides with the order 
parameter amplitude of the two-dimensional attractive 
Hubbard model at its maximum T c (as obtained using 
MCMF). In other words, for large U the bilayer achieves 
its optimal T c by adjusting the inter-layer tunneling t± 
to a point that maps the bilayer onto a single layer at- 
tractive Hubbard with an optimal U /t ratio. Since the 
average pairing amplitude of the latter is of order t Eq. 
(JM| implies the ij_ imax ~ VtU scaling mentioned above. 

Finally, in order to obtain a rough idea of the effect of 
disorder on T c . max we resort to the mean-field treatment 
of Ref. [IJ. We apply it to a bilayer with U = t where 
in the absence of disorder it differs from the QMC and 
MCMF results by a factor of 3, see Fig. H The mean- 
field approximation consists of decoupling the interaction 
term -UfLftff^fn -> -A|/ if / 4 - Aj/^/^, and solv- 
ing the BdG equations for the disordered bilayer at finite 
temperature. This is done under the self-consistent con- 
dition Aj = U(fi^fi^) and with a phase twist in the x 
direction, which enters the kinetic part of the Hamilto- 
nian according to 

-* E £ e*C*<-*')/ a °cUv. (37) 

The free energy, F, calculated from the BdG solutions, 
is then used to evaluate the bare phase stiffness p° s (T) = 
(a/L) 2 d 2 F/dq 2 , which includes the physics of thermally 
excited quasiparticles but ignores the renormalization of 
the stiffness by phase fluctuations. 

Figure [6] depicts the estimated maximal critical tem- 
perature obtained from the approximate BKT criterion 
T c ,max = (7r/2)p°(T Cimax ,£_ L ,max) for a bilayer with L = 
10a, and n — 1.5. The results were averaged over up 
to 200 realizations of disorder in which every Vi was 
drawn independently from a uniform distribution in the 
range [— V, V}. As shown, the disorder has little effect on 
7c !ma x as long as it is smaller than the metallic hopping 
amplitude. For V/t > 2.5 we found that the system's 
dimensionless conductance, g, obeys g < 1, indicating 
that strong localization effects become important in this 
range. Such effects enhance phase fluctuations and are 
expected to induce a transition to an insulating phase, 
which is not seen in the small system considered by us. 
We also found that the value of f_i_, max is not affected by 
the disorder and remains fixed at i_L !max ~ 0.48t. 



0.1 




1 2 3 4 5 6 

V/t 

FIG. 6: The estimated maximal critical temperature T c , max 
deduced from the approximated BKT criterion T Cimax = 
(•7r/2)p2(7c,max) as function of the disorder strength V for 
U/t = 1 and n = 1.5. The bars depict the standard deviation 
of the distribution of T c . max . Based on our results for the 
clean system we expect the actual maximal T c to be about 
3 times smaller than the values presented here. For V/t > 2 
strong localization effects likely cause T c to drop even faster 
and to eventually vanish at a critical disorder strength. 

IV. CONCLUSION 

Superconductivity contains within itself a built-in ten- 
sion between its two essential prerequisites: pairing and 
phase coherence. In all known examples where interac- 
tions are strong such that they lead to a large pairing 
scale a concomitant suppression of the superfluid stiff- 
ness, and with it of T c , takes place. Canonical exam- 
ples are the strong interaction regime of the attractive 
Hubbard model and the underdoped cuprate supercon- 
ductors. It was suggested that a possible way out of 
the dilemma is to separate the pairing medium from the 
conduction electrons, thus maintaining their large phase 
stiffness^ This strategy was pursued in Ref. LU and its 
qualitative feasibility was demonstrated. However, the 
consequences of superconducting phase fluctuations on 
the induced stiffness at the interface between the two 
subsystems were not considered. 

In the present work we have shown that phase fluctua- 
tions lead to an exponential suppression of the superfluid 
stiffness when the interlayer coupling is small. For inter- 
mediate coupling these fluctuations cause a reduction of 
T c by factor 3 relative to its value based solely on ther- 
mal excitation of quasiparticles. We believe that our re- 
sults have broader implications beyond the specific model 
studied here, as they point to the fact that renormaliza- 
tion of the stiffness by phase fluctuations is important in 
systems with long range phase couplings. 

Despite the renormalization, long phase couplings are 
preferable to short ones from the perspective of increas- 
ing T c at the surface of a phase fluctuating superconduc- 
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tor. As we showed, when a metal is overgrown on the 
surface of a superconductor, additional phase couplings 
between sites on the interface are established up to dis- 
tances of the order of the metallic thermal length. This 
is true whether the added metal forms a two-dimensional 
layer or is thicker than the thermal length and can there- 
fore be considered as three dimensional. However, within 
this range the decay of the coupling, technically described 
by the cooperon, follows r~ d in the (i-dimensional case. 
This is a result of the fact that in three dimensions the 
electrons that mediate the coupling spend part of their 
time moving in the direction perpendicular to the sur- 
face, thereby lowering their probability to reach longer 
distances along the surface within the allotted coherence 
time 1/T. Consequently, our analysis suggests that one 
should attempt to make the metallic layer as thin as pos- 
sible in order to maximize the induced phase couplings 
and with them T c . 

Our results show that disorder has little effect on T c , 
as long as it is weak enough not to cause strong localiza- 
tion. From a practical point of view, a particular form of 
disorder, namely interface roughness, may actually ben- 
efit the enhancement of T c £2- This statement stems from 
the fact that when the pairing scale is of the order of 
the metallic bandwidth maximal enhancement occurs at 
values of the interlayer tunneling which are comparable 
to the metallic hopping amplitude. Such strong tunnel- 
ing across the interface may be difficult to achieve. An 
example can be found in the cuprate bilayerss^ where 
the intrinsic in-plane hopping is much stronger than the 
inter-plane one. Nevertheless, if the interface is not per- 
fect so that the pairing medium and the metal locally 
interpenetrate each other electrons may tunnel between 
the two laterally, exploiting the large hopping amplitude 
in this direction. 



where 



IJ u UJ m ,U) n j 



Vy 0) (w m ,u;„) 

r(l)/ 



P 



V^\u) m ,u) n ) = -t 2 x ^2 [(I + <T3)Gij(uj m - Q n )Fi(Q n ) 
xF*(tt n + uj u - w m ) - (I - a- 3 )Qji(-u} n - Sl^Fjtrin) 



(Ul m — UJ„)8i(uj m — LO n )(T 3 SijS u 



xF*(tt n + LJ n - ui m )] 



(A2) 



Here, and throughout, w„ and fi„ are fermionic and 
bosonic Matsubara frequencies, respectively. We also de- 
noted 



Fi(r) 



,-i6i(r)/2 



^Fi(n n )e 



-iSl„ 



(A3) 



Our goal is to integrate out the fermions to obtain 
S e = -Trln[Af(°) + V W + V {1) ] 

« -Trln[M(°)] - Tr {m^~\v^ + V™]} 

+ iTr{Af(°)~V 0) +^ (1) ]} 2 • (A4) 

A/' ) is independent of 9 and thus plays no role here. It 
is also trivial to see that the linear term in vanishes. 
Hence, we turn our attention to 



-Tr 



8t \ t 2 ± N F a 2 



min(l,— 



(A5) 



evaluated for fi„ < Ao. This, however, is negligible com- 
pared to 
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Appendix A: Derivation of the phase action 

Here we provide details concerning the derivation of 
the phase action, Eq. As discussed in Section Hi Bl 

we assume that the superconducting amplitude is fixed 
at its zero temperature value |Ai(r)| = Ao and pro- 
ceed to write Eqs. (|2I3|) in frequency space. In terms 
of *J(w m ) = [// t (w m ), fn(-LJ m )] they become 



(Al) 



iTr U./ (0r V (0) M(°)~V (0) 



1 Y,^m^\ 2 +°(— 1 



8TA 



•fi 



(A6) 



which constitutes the first term in Eq. 

The term |Tr [m^'V^M^T V (1 > 1 contains two 
contributions. The first is 



Ci 



E E E 



t 4 A 2 



x Fi(Q n )Fj (fl n +u m ~ u n )Gij{-u v 

in 



(^+^ + A 2)(^+ M 2 + A 2 ) 

•fin) 

- -m 



(A7) 



Its evaluation requires the disorder average Pij (u) n , fi„) = 
(Gij (w„) Gij (Wn + fin)}- The metallic layer is assumed 
to contain Gaussian disorder obeying (V*) = and 
(YiVj) — nV 2 8ij, where n is the average number of im- 
purities per site and V 2 = {V 2 )/n. The disorder scatter- 
ing is characterized by the elastic mean free time r e = 
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[27rA^a 2 nV 2 ] -1 and is assumed weak, kpl e ^> 1. The 
leading contribution to the above average comes from the 
cooperon, i.e. , the sum of ladder diagrams whose legs 
are constructed from disorder-averaged Green's functions 
and whose rungs are disorder lines. In momentum space 
the cooperon takes the form= 



1- (2^a 2 r e )- 1 P (Q,w„,«„)' 

(A8) 



Here 



1 



-Po(Q, Win fin) — T7 



1 



1 



+ fi„) - £ k +Q + 2^sign(w„ + n n ) ' 



(A9) 



with £k = — 2i[cos(fc x a) + cos(fcj,a)] + e — /i, reduces to 
the cooperon of the clean system in the limit r e — > oo. 
For Q = |Q| -C kp it becomes 

o/n n A >u 2 |sign(tj» + »») -sign(o; n )| 

V(l + |fin|r e ) 2 + (QZ e )2 

(A10) 

For a diffusive system, and as long as Ql e <C 1 and 
|w ra |T e <C 1, one may expand the square root in Eq. 
(|A10p . and plug the result into Eq. (fAlj) to find 



2 |Sign(w„ + fi n ) - sign(w ri 



(All) 

where D = l 2 /(2r e ) is the diffusion constant. At this 
point we would like to note that an additional time scale, 
the phase breaking time r„, can be phenomenologically 
introduced into the problem and act as a mass term for 
the cooperon. It is knowr*2£ that interactions in two di- 
mensions lead to t^ 1 ~ (T In g)/g, where g is the di- 
mensionless conductance of the disordered layer. Our 
treatment of the disorder is valid for j> 1 and therefore 
t~ -C T. Consequently, it does not alter our results and 
is not considered further. 

Transforming back to real space leads in the clean limit 
and ru ^> fcp 1 to 



Pij(u> n , fin) = -^Npa 4 |sign(w„ + Q n ) - sign(w„)| 



1 



-\fim\rij/VF 



(A12) 



VFTij 

while for the diffusive system and 3> l e the result is 
N F a 4 



2D 



x K [r 



|sign(w„ + Vl n ) - sign(w„)| 




(A13) 



where Kq(x) is the Bessel function of the first kind. 



Plugging Eq. (|A12[) into Eq. ([ATI) , defining u = 
fim — w m , uj' = f2„ + u) m and carrying out the summa- 
tion over uj m and uj n using the fact that at the relevant 
low temperature range they can be approximated by in- 
tegrals one obtains C\ = J2i j where in the clean 
limit 



t\N F a A T 



4^2 r P 4 



n dr - 



e -A |ri-T2| e -Ao|T3-r 4 | 



OVFTij Jo " 

x F^Fiir^i^F^n) |sign(w') - sign(w)| 



x e 



i[w(Tl — T 3 )+a)'(T 2 — Ti)] — U\rij/VF 



(A14) 



For Tij > lp = vp/T the sum is dominated by the lowest 
frequency difference and decays as (l/r i j)e~ 2 '* rij / lT . 
For nj <C It the sums may be approximated by integrals 
with the result 



n _ t 4 ± N F a 4 



8ir 2 VFrij 
(n- 



4 

n ^ 

n=l 
T3)(r 2 



e -A |ri-r 2 | e -Ao|r3-T4| 



-4) + (^) 2 



[(r 1 -r 3 ) 2 + (^) 2 ][(r 2 -r 4 )2 + (^)2] 
x F^F^F^t^F*^). (A15) 



Owing to Eq. dXBJ) (Fi(T)F?(T')) 



-A |t-t'|/2 

) 

which means that may be considered constant within 
a time window of order Aq 1 . This fact and the expo- 
nential factors in Eq. (|A15|) enable us to approximate 
Fi(n) = F,(r 2 ) and F*(t 3 ) = F/(t 4 ), leading to the 
coupling cos[^i(r) — 0j(r')] mediated by pair tunneling 
between the two sites, and constituting the second term 
in the action, Eq. (0}. The associated kernel may be 
evaluated in the range a <C r.y <C £ = vf/Aq by taking 
rij — > in Eq. ()A15|) and carrying out the integration 
over t 2 and r 3 . The result 

,4 



X(r,r) = 



i 4 A^a 4 
8n 2 vpr 



° r Ei[— Aqt] - e~ Aor Ei[A r])' 



,A r 



+ 7Te 



2„-2A |t| 



(A16) 



where Ei(a;) is the exponential integral function, is ap- 
proximated by Eq. ©. For lp ^> > ^ we may 
expand the integrand in Eq. (IA15[) to second order in 
Ti — t 2 and T 3 — T4. Integration over t 2 and T3 then yields 
Eq. ©. 

Repeating the derivation using Eq. (| A13|) for the dif- 
fusive case produces a similar coupling with the kernel 



K{r,r) 



t 4 ± A N F a 4 T 



1 



u)>0 v 



A 2 



Kr 



x r^ e - 2A °i"i sin[2a;(r 

./-oo T + V 

x [ui cos(2wr7) + Aq sin(2w|?7| 



■TI)] 



(A17) 



The exponential decay of K$(x) for x > 1 induces the 
decay of the kernel when > lp = y/D/T. This 
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also means that for It ^> r,j ^> £ = ^/D/Aq only 
w < Ao(£/r) 2 contribute. Since the exponential factor 
in Eq. (|A17|) implies \t)\ < Ao we may expand the in- 
tegrand in a; 1 77 1 < 1 and evaluate the integrals with the 
approximated result Eq. (0). 

The second contribution to the effective action coming 

from ±Tr Uf(°)~ V«M(°)~ V«~ 



Co = 



E 



M 2 + Ag)(o;a+ / i 2 + Ag) 
x Fi(Q n )Fj(Cln + u m - ui n )Qij{-ui m - Q n ) 

X Fj(Cl m )F* (fi m +U„ - U m )Gji(-U)n - Hm)< (A18) 

Following a similar line of derivation to the one taken 
above we find that this contribution corresponds to a 
coupling cos {| [#i(/r) — 0j(t) — ^(t') + 9j( T ')]}, induced 



by single-particle tunneling between the sites£L However, 
we find that it has a negligible effect in comparison to the 
pair-tunneling term. For example in the clean system and 
for It Tij £ the associated kernel reads 



- t*JV F oV 
K{r, t) = 



1 



- {r/vF? 



27t 2 A 4 uf r [t 2 + (r/v F ) 2 



(A19) 



Applying the mean field treatment of Sec. Ill Dl to 
this term yields an average coupling constant g ~ 
(2t 4 L N F a 2 ^ 2 )/(A r rod 7rA^) to be compared with g — 
(tiN F a 2 /N rod A 3 Q )\n(lT/0, Eq. 



Finally, the Tr 



piece in Eq. 



(|A4j) can be shown to result in a (d T 9i) 2 term, which is 
of order t\ and is therefore insignificant compared to Eq. 
[1. 
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